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We study discrepancy between the analytical definition and the numerical implementation of the 
McLerran-Venugopalan (MV) model. The infinitesimal extent of a fast-moving nucleus should retain 
longitudinal randomness in the color source distribution even when the longitudinal extent approx- 
imates zero due to the Lorentz contraction, which is properly taken into account in the analytical 
treatment. We point out that the longitudinal randomness is lost in numerical simulations because 
of lack of the path-ordering of the Wilson line along the longitudinal direction. We quantitatively 
investigate how much the results with and without longitudinal randomness difi^er from each other. 
We finally mention that the discrepancy could be absorbed in a choice of the model parameter in 
' the physical unit, and nevertheless, it is important for a full theory approach. 

^ \ PACS numbers: 24.85.-|-p, 12.38.-t, 25.75.-q 
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^ \ I. INTRODUCTION 

It is a tempting idea to approximate the wavefunction of an interacting particle by mediating gauge fields sur- 
^vq rounding the particle. Such a description of a fast-moving charged particle by virtual photons is long known as the 
Weizsaker- Williams approximation. If we apply this idea to the strong interaction, we can approximately give the 
wavefunction of a heavy hadron by gluon fields which extend from the local color distribution inside the hadron. 
The formalism of the Color Glass Condensate (CGC) can accommodate the non-Abelian extension of the Weizsaker- 
Williams approximation [l|, Q • 

The CGC formahsm aims to embody the parton saturation picture arising as a result of the small-x evolution at 
^ high energies [Sl] . The geometric scaling beautifully indicates the existence of the saturation scale as a function of 
Bjorken's x The CGC formalism has been forming an important building-block of Quantum Chromodynamics 
(QCD) at high energies. With a Gaussian approximation for distribution of the random color source and its dispersion 
given by the saturation scale [H], we reach the McLerran-Venugopalan (MV) model [1]. 

So far, the MV model has been well examined analytically in the case of a light projectile (e.g. a color dipole, a 
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, i.e. high-energy nucleus-nucleus collision [1 
^vq ■ analysis is limited to the numerical method 
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proton, etc) scattering off a dense CGC target 
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conditions) are analytically known in a simple form 



[Id . lUl, llJl- The scattering problem of two CGC objects, 
|17| is hard to solve analytically, however, and the possible 
I22I . Only the initial fields right after the collision (i.e. initial 
llSl . [23} . The initial chromo-electric and chromo-magnetic 



I . fields are thus calculable by means of the proper Gaussian average over the color source distribution [17|, l21| . 
J — . In this work we will address a problem that the numerical formulation of the MV model assumes a crude ap- 

, proximation which leads to a significant difference from the right answer. We will clarify the physical meaning to 
ILJ ' ascertain that the assumed approximation is not acceptable. The crucial point is that the formulation involves two 
distinct limits in the longitudinal direction; one is the vanishing extent of a nucleus and the other is the vanishing 
correlation length in the random color distribution. We will then quantify how much the numerical procedure turns 
\ out to underestimate the expectation value of the initial fields as compared to the analytical answer. 

One might wonder if the issue to be discussed here is academic and only a technical detail, for phenomenological 
models are useful as long as they can fit in with empirical data. The existing numerical simulations in the MV model 
are, in fact, fairly consistent with experimental observations at RHIC (Relativistic Heavy-Ion Collider). We would 
emphasize, however, that the MV model is not a fitting model but a theoretical description based on QCD. Therefore, 
one should not make little of correctness checking as a theory problem. 

Even though the numerical calculations underestimate the field expectation value, one might suspect that its effect 
is negligible from the fact that the numerical simulation agrees with RHIC experimental data. We shall see, however, 
that the approximation used in the numerical simulation causes a factor over 10 in the initial energy density. Such 
a discrepancy is not visible in the former numerical works in the MV model because the choice of the MV model 
parameter can absorb it. All the observables scale in accord with /i in the saturation regime, and in the numerical 
simulation usually, /i is determined so as to reproduce the particle multiplicity. As a result, for instance, twice larger 
/i leads to 16 times larger energy density which may compensate for the discrepancy in the dimensionless coefficient. 
This "posteriori reasoning" would sound fine as a phenomenological approach. The MV model is a theory, however, as 
we emphasized above. The essential notion of the CGC idea is that a universal description for the hadron wavefunction 
becomes possible once the saturation occurs. Hence, /.t could be fixed independently by information obtained in the 
pQCD calculation or in the deep inelastic scattering (DIS) through the relation to the saturation scale Qs as a function 
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of fjL. [See a recent nice work [IJ] for careful discussions.] It is, of course, quite difficult to specify for RHIC in this 
way because the relevant Qg at RHIC is not unique from the kinematics unlike the DIS case. 

The reason we would stick to rigorousness also lies in the time scale of early evolution after the heavy-ion collision. 
A larger /i may lift the underestimated energy and multiplicity up at the cost of changing the time scale oc 
Moreover, our finding must have a significant effect on the "Glasma" instability proposed in Ref. [Hf. Because the 
correct treatment in the longitudinal direction is important in our argument, as we will see later, the instability with 
respect to the longitudinal variable (rapidity) should naturally be affected. Actually the instability is found too weak 
in Ref. [1^. It is likely that the numerical simulation in Ref. [1^ underestimates the instability strength in the same 
way as for the energy density. If so, there might be a chance that the genuine Glasma instability would occur faster 
to contribute more or less to the early thermalization mechanism. 

One more point that motivates us to reveal the discrepancy stemming from the numerical approximation is the 
infrared (IR) property of the MV model. The IR cutoff has a physical meaning, while the ultraviolet (UV) cutoff 
should be zero in the end. The dependence on the IR cutoff strongly relies on the approximation. This feature would 
bring about uncertainty, as is discussed in Ref. [131 > in addition to the overall dimensionless factor. 

In this work we will adopt the same prescription for the analytical calculations as in the numerical simulation to 
regularize the UV and IR singularities; we use the lattice formulation with the spacing and the system size given by 
a and La, respectively. This enables us to make a direct comparison between the analytical and numerical results. 



II. MCLERRAN-VENUGOPALAN MODEL 



Here we will make a quick review on the MV model applied for the relativistic heavy-ion collision (two-source 
problem). If we have two sources which represent particles moving fast in the positive and negative z-directions 
("right" and "left" respectively), we can write the corresponding current using the light-cone coordinates as 



(1) 



where we dropped x+-dependencc in the right-moving source and a;~-dependence in the left-moving source as usual 
because of the Lorentz time dilatation. Then, converting the coordinates from the light-cone variables to the 
proper-time t and the rapidity rj which are related by x^ = {t / \/2)e^^ , we can express the initial fields right after 
two nuclei collide as [13, Sal 



Un - ,-„rr^(l) ^(2)1 r (2) (1)1 



(2) 



and the transverse fields (i?*, B') are zero [2l|, under the assumption that '{x±,x ) and ' {x±, x'^) behave like 
S{x~) and S{x~^) respectively near the light cone This assumption is the case because of the Lorentz contraction 

(n) 

when two nuclei travel at the speed of light. In the above, the color matrices, s, are the gauge field configurations 
in the radial gauge {Ar oc x^A'^ + x^A^ ~ 0) created from each color source p*^"'. The solution must satisfy 



(") 



with a constraint that takes a pure gauge form so that its field strength is vanishing. 



The covariant gauge is most convenient to get an explicit solution. By rotating the solution obtained in the covariant 
gauge to the radial gauge, we can write down the desirable solution as {§]; 



a. 



i'\x^,x~) = --V{x^,x-)d,VHx^,x'), a'f\x^,x+) = --W{x^,x+)d,W\x^,x+) 



(2), 



(3) 



where the gauge rotation matrices are 

) = T^x- exp 
W\x^,x'^) = Vx+ exp 



^9 



«5 



/ dz I d^y^Ga{x^-y^^p^^\y^,z ) 
f dz+ [ d^y^ Goix^-y^) ~p^^\y^. z+) 

J —00 J 



(4) 



Here V^i denotes the path-ordering with respect to x^ . We remark that p in Eq. ((4]) is the color source in the 
covariant gauge that is not identical to original p in the radial gauge. Since the Gaussian weight is a gauge-invariant 
function of tr[p2] = tT[p^] (see Eq. ([7])), we do not have to discriminate p and p practically. We will thus not use 
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the tilde any more but just write p regardless of gauge. The two-dimensional massless propagator is defined by 
c?5_Go(a;^) ~ — J*^^^ which is singular in both the IR and UV sectors. Let a and L be the lattice spacing and the 
number of the lattice sites, and then we can express G'o(£Cj^) in the lattice regularization as 



L/2 

Eexp 
^ 2 — cos(27mi/L) — cos(27rrt2/i) 



[i{xil'Kni/ La + X22.'!in2 / La)] 



(5) 



The summation is supposed to exclude a singular point rti = n2 = because we impose the global color neutrality 
{p^ = 0) = 0, so that the singular point ni = n2 = does not appear at all in Eq. (jH). It is easy to check that 



(n) 



■P^/[(27r) p]^], in the contimmm limit where 



P 

Eq. ^ is reduced to the standard expression, Go(a;j 
a — > and L — s- cx) are taken. 

The later-time dynamics is uniquely determined by the equation of motion with the initial fields given by Eq. 
Any physical observables are therefore given in terms of V and W in principle, which we shall denote as 0[V,W] 
generally. We can compute the expectation value by taking the average; 



{0[V,W]) = / [V'][V^]WW[pW]>V(2)[p(2)]o[y^H^]] 



(6) 



with the Gaussian weight, 



W(2)[p(2) 



exp 



exp 



dx (fx 



tr[p(i)| 



a;_L, X 



dx^(Px± 



tr[p(2)(a 



(7) 



The normalization of the color trace is understood as tr[i™i"] = -^(5™", where t™'s are the SU(iVc) algebra in the 
fundamental representation. The scale /l(^^'(2;~) and thus p*-^^ (a;j_, x^) (or /i(2)(x+) and thus p^'^\x_i,x~^)) may have 
finite extent in the x~ (or x~^ respectively) direction since the Lorentz 7-factor is finite in fact and also because the 
small-x evolution by the JIMWLK equation should spread longitudinally in the color distribution. 

If we are interested in evaluating the initial energy density in the heavy-ion collision, 0[V, W] should be tr[(£''')2 + 
that is [i3,[2l, 



Co 



(8) 



where we introduced a notation to indicate the diagonal component, i.e. (afa'j) — SijS°'''{aa) for a*^^-* and a'-^^ 
respectively. We will make use of this notation again when we show the numerical results later. 



III. QUESTION - NUMERICAL APPROXIMATION 

In the numerical implementation, for practical reasons, it is difficult to compute the expectation value of the 
Wilson lines ^ as they are. In the first approximation the longitudinal finiteness is only negligible as compared to 
the transverse size and one can take the following limit; 



(X^ ,X-)^ (X^) S{X~ ) , p(2) (^^ ^ ^ -(2) ) ^(^+) 

Then one might anticipate that all physical quantities are given as a function of the integrated scale, 

[mW]'=/ dx-[p('Hx-)]\ [fl^'^]'= dx+[^^('Hx+)]\ 



(9) 



(10) 



As a matter of fact, we can verify this expectation in the analytical evaluation associated with only one source [8l[9l.[l^. 
that is, only [/i^^-']^ (or [/i(2)]2) appears in the Gaussian average of a function of V (or W) with the weight yy^^^[p^^^] 
(or >V(2)[p(2)j respectively). 
In the popular numerical formulation the Wilson lines ([4]) simplify approximately under the limit Q as 



V^x^,x-) V^{x^,x-) = exp 
W''{xjL,x+) ^ W\x^,x+) = exp 



ig Jd^y^ Go{x^^y^)p^^Hy^)0{x-) 
ig [d^y^Goix^-yJp^^\yjeix+: 



(11) 
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These would be, of course, exact if there were not for the path-ordering Vx± or if in the Abelian gauge theory. 
The question we are addressing in the present paper is the vahdity of this naive prescription (jlip that is commonly 
employed in the numerical simulation. 



IV. MATHEMATICAL REPRESENTATION 



Let us focus only on the ^-sector in what foUows because the initial energy density is factorized into the t^-sector 
and the H^-sector as is manifest in Eq. (|5]). Needless to say, exactly the same argument should work for the M^-sector 
as well. For notational simplicity, then, we shall omit the superscript (1) which represents the right mover. 

The subtle point comes from the fact that another Dirac delta function is involved implicitly in the formalism 
besides Eq. ([9]). Namely, the Gaussian weight ([7]) leads to the following correlation function; 

{pa{x^,x-)p^,{v^,y-)) = {ii(x-)\ ' J"" b{x^-v^) 5{x- -y-) . (12) 

We need to deal with two 5{x~y& properly in order to formulate the problem in question in an appropriate way. For 
this goal it is convenient to introduce a regularization to the longitudinal Dirac delta functions. We modify the Wilson 
line accordingly as 



V^{x^,x ) = V^- exp 



J — OC J 



(13) 



where the regularized color source must be a smooth function satisfying 



\\jnp^{x^,x )=p{x±)5{x ). (14) 
We shall introduce another regularization for the correlation function in a way as 

{paix^,x-) pb{y^,y^))^ = [pix^)] ^ S'''' S{x^-y^) d(;{x- -y-) , (15) 
such that the regularized delta function must satisfy 

\iiR6(;{z-) = 5{z-) . (16) 

We are now ready to elaborate the question in a mathematically sophisticated manner. The replacement of Eq. (jlip 
is justified if e comes first before C, ^ 0. In the analytical calculation, on the other hand, the relevant limit is 
C — > followed by e — > later. Thus, a mathematically sensible description of the question pT|) should be 

lim lim/OfFj)^ = \im\im(0\V,]),. (17) 

Here, the left-hand side corresponds to the numerical irnplementation pT|) and the right-hand side corresponds to the 
one that is commonly assumed in the analytical works 0, [^, [lol | . 

We sketch the intuitive interpretation of e and C in Fig. [TJ Roughly speaking, e is the longitudinal extent of a 
fast-moving nucleus and C is the correlation length of the color distribution inside the nucleus. The physical limit 
should keep e > C as is the case in the right-hand side of the question (fT7)) . If e goes to zero first, the longitudinal 
structure of randomness is lost. Then only one infinitesimal "sheet" of the two-dimensional random color distribution 
is left and the path-ordering becomes irrelevant (see also Fig.d]). The numerical prescription (jlip hence implies such 
unphysical ordering of two noncommutativc limits. 



V. COMPARISON 



We will explicitly confirm that the order of two limits is noncommutativc indeed. We will begin with the simplest 
example of the tadpole expectation value. Then we will proceed to the case of the gauge fields which is directly related 
to the estimate for the initial energy density by Eq. ([8]). We note that the tadpole gives the scattering amplitude 
between a quark in the color fundamental representation and a CGC medium. This scattering amplitude is IR screened 
by long-ranged color interactions. In contrast to that, the gauge field correlation function involving four Wilson lines 
contains a color-singlet component which is free from screening. However, in this case, the derivative acting onto the 



FIG. 1: Schematic picture of a nucleus with two regulators e and Intuitively e signifies the longitudinal extent of the whole 
nucleus and stands for the longitudinal correlation length inside the nucleus. If randomness of the color source distribution 
is attributed to color confinement in each nucleon, corresponds to the longitudinal extent of the nucleon. 



Wilson lines is UV harmful leading to singular behavior. It is possible to take the limit of vanishing UV cutoff (a ^ 0) 
at finite time (r ^ 0) owing to non-linear evolution as has been closely investigated numerically in Refs. [2ll . [2^ and 
analytically in Ref. [I4I. The initial energy density in the heavy-ion collision at t = is, in this sense, an ill-defined 
quantity in the a —^ limit. We will estimate it here, nevertheless, to compare the analytical and numerical outputs. 
It is partly because, whether it is a physical quantity or not, we are pursuing to exemplify a significant discrepancy 
anyhow. We will do this also partly because we know from Ref. that the logarithmic singularity ~ [ln(La/a)]^ 
at T = translates into ~ [ln(La/r)]^ at r ^ with the overall coefficient unchanged. Thus, if we find a difference 
in the overall factor of the initial energy density as we will do below, the UV safe energy at r 7^ should receive 
the same overall factor. Accordingly the initial energy estimate here could serve as correctness checking for UV safe 
observables indirectly. 



A. Tadpole 



In the simplest case of the tadpole operator, i.e. 0[V] = V\ we can perform a quick comparison even without 
resorting to the numerical method. We already know the analytical answer for the right-hand side of Eq. (fTT]) . That 
is given by 



lim lim(V',^)^ = exp 



-9 M 



L(0,0) 



where L(0,0) is the notation in Ref. ^] which is defined as 

2 L/2 



L(0,0) = J (fx^Go{x^)Go{xj 



a 
4i2 



E 



n; = l-L/2 



[2 — cos(27rni /L) — cos(27rn2 / L)] ' 



2tt J 



2tt 



(18) 



(19) 



in the lattice regularization. Here, as in Eq. ((5]), the zero- mode (ni = n2 = 0) is to be removed by neutrality. The 
quadratic form approximates the sum quite well with a coefficient 0.962 that we find numerically. 
As for the left-hand side of the question (fTT)) . we have to perform the following Gaussian integral; 



lim lim(V"/) 



C 



[dp] exp 



«5 d'^y±Go{xj_-y^) p{y^) 



exp 



d^x 



tT[p{xj 



9 fJ' 



(20) 



It should be mentioned that the first exponential part is a matrix, which makes the Gaussian integral hard to 
accomplish. Although it is a tough calculation for arbitrary SU(iVc) group, the SU(2) case {Nc = 2) is feasible 
immediately because the SU(2) exponential matrix is easily manipulated. After some calculations we find that the 
above Gaussian integral results in 



lim lim(y/). 



1-5^7^(0,0) lexp 



which obviously differs from Eq. p8|l with Nc = 2 substituted 

lim lim(V/(a; 1 = exp 



-9Vgi(0,0) 



(21) 



(22) 
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FIG. 2: Tadpole expectation value as a function of L with the parameter choice g^Jia = 0.17. The left figure is for the SU(2) 
group and the right for the SU(3) group. The solid curves represent the analytical results by Eq. (|18|l . The open circles 
represent the numerical results by the Monte-Carlo integration. The thin dashed curves show Eq. (|21|l in the SU(2) case and 
Eq. (|25p in the SU(3) case, respectively. The SU(2) and SU(3) numerical data agree well with the thin dashed curves, which 
means that our Monte-Carlo integration works nicely. 



It is interesting to see from the above that Eq. (j21|l can be a good approximation to Eq. (j22|) as long as g jlaL is 
small enough to allow for the Taylor expansion. 

It should be instructive to compute Eq. (|20p numerically in the Monte-Carlo integration though we already have 
the answer. In order to compare to the existing numerical calculations in literatures, we shall adopt the parameter 
choice same as used in Ref. [22l|. That is, 



9 1 2- ni7 
— = - , g ^ia = 0.17, 

47r TT 

was chosen in Ref. [23] and then the nucleus size is given by 



L = 700 . 



(23) 



Ra = 



aL 



(24) 



We change the value of L to see the functional form. Therefore, if we increase L while keeping g^fia fixed at 0.17, the 
nucleus size or the infrared cutoff Ra grows up in proportion to L. 

We plot the analytical formula (fT8|) as a function of L by the solid curve for the SU(2) case in the left of Fig. [2] 
and the SU(3) case in the right of Fig. O Since g^fla is a constant, what Fig. [2] means is the i^^-dcpendence of the 
tadpole expectation value. The open circles in Fig. [2] represent the numerical results by means of the Monte-Carlo 
integration. We took 200 ensembles to calculate the expectation value. In the SU(2) case the numerical results agree 
well with the expression (|2ip . Also, an analytical expression. 



lim hm/yt). 



C-*Oe-^0 



1 - .9Vii(0, 0) -I- gV_Li(0^ 0)2 ) exp 



gVgi(0,0) 



(25) 



can nicely fit the SU(3) results. In any case, it is obvious that the numerical results deviate from the analytical 
formula (fT5|) substantially. 



B. Gauge Fields 

In view of the tadpole results, one might have thought that the discrepancy is only minor. The deviation may look 
small, however, simply because the expectation value of color non-singlet operators is exponentially suppressed by 
the system size Ra- This fact becomes manifest once we consider some other operators that contain a color singlet 
component. 

Here, let us elucidate a more complicated situation than the tadpole, that is, the expectation value of gauge fields 
(aa) (see Eq. (|S]) for our notation) which have a contribution from the color singlet. In this case the Gaussian integral 
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FIG. 3: Gauge field expectation value as a function of L. The analytical results are given by Eqs. ((26} and (|27|l . The parameter 
g^fia is fixed to be 0.17 for the a-fixed results and g^fia scales as 0.17 x 700/1/ for the -RA-fixed results. The left figure is for 
the SU(2) group and the right for the SU(3) group. 



is too tedious to accomplish, so we will rely on the numerical Monte-Carlo integration only. On the other hand, the 
analytical calculation with the path-ordering is still possible and it follows; 



(aa) = l(5'/2)2a, (26) 

r 



where 



1 1 1 1 



ln(1.36L) = - ln(2.41i?A/a) , (27) 



4L2 2 - cos(27mi/L) - cos(27r7i2/L) 47r ^ ' ' A 

which does not depend on Nc at all and monotonically grows up with increasing L. We numerically find a constant 
1.36 in the logarithm which approximates the sum. This expression is, hence, both IR and UV singular. That is, 
{aa) — !■ oo whenever L oo (i.e. either Ra — » oo with a fixed or a — > with Ra fixed, see Eq. (|24p ). More 
importantly, the analytical formulae (pS]) and claim that the IR behavior as Ra oo and the UV behavior as 
a — > are completely identical. This property is, however, no longer the case in the numerical results with wrong 
approximation assumed, which could be a possible explanation for IR stability found in Ref. [l8[. 

Figure [3] shows the results by the Monte-Carlo integration. The shape of {aa) / [g"^ jxf as a function of L clearly 
depends on whether Ra oc L increases with a fixed or a cx 1/L decreases with Ra fixed. In the case that a is fixed, the 
calculation goes just in the same way as in the previous subsection; we choose g'^fia = 0.17. When we keep Ra fixed, 
we adjust the lattice spacing as g^fla = 0.17 x 700/L, so that g'^Jia becomes 0.17 for L = 700, which is completely 
equivalent to the procedure in Ref. [22]. Thus, the a-fixed curve and i?^-fixed curve should meet ai L — 700, as is 
certainly the case in Fig. [31 

The numerical results hardly depend on N^, the left and right figures of Fig. [3] look almost the same, though the 
SU(3) results are shghtly smaller than the SU(2) ones. The numerical calculation leads to g"^ (aa) / {g"^ fi)'^ — 0.194 at 
L = 700 for the SU(2) group, while the SU(3) group results in 0.150 at L = 700. Differently from the tadpole case, 
the numerical data underestimate the expectation value that is g'^ (aa) / (g'^ jl)'^ ^ a — 0.546 in the analytical method. 
In the next section, we will discuss how we may be able to remedy this situation. 



VI. IMPROVEMENT 



We can improve the situation by inserting the infinitesimal sheet in longitudinal extent in a way as sketched in 
Fig. m We denote the number of sheets by N^j. We note that is not the longitudinal coordinate as in Ref. but 
it is the number of slices within infinitesimal extent. Hence, all the numerical results we have seen so far correspond 
to N,j — 1. We can recover the full path-ordering in the Njj oo limit. 
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FIG. 4: Schematic picture of how to improve the numeri- FIG. 5: Tadpole expectation value for the SU(3) case at 

cal results. Each blob represents the two-dimensional sheet A^,, = 1 (same as shown previously), A^^ = 2, and A'^^ = 3. 

without longitudinal randomness. In the A^^ oo limit 
we can recover the full randomness structure in longitudinal 
extent which is infinitesimal in the e — > limit. 



For example, in the tadpole case, we shall modify the Gaussian integral in Eq. ([20)1 into the form of 



(^')^ =n h^pn] exp 



exp 



- (fx 



(28) 



to retrieve the analytical results in the limit of iV^,, 
Eq. PH)) . It is, however, impossible to take the Njj - 



oo. We have to divide ff by iV,, 
c» limit in the practical procedure. 



to make it consistent with 
Instead, in this section let 



us focus only on N^j = 2, iV^ — 3, and — 10 to demonstrate the tendency of how the numerical outputs could move 
from the structureless N.^ = 1 results toward the analytical answer. We did calculate in the SU(2) case as well as in 
the SU(3) case, but we will present only the SU(3) results here, for the gauge group makes no qualitative difference. 

The improvement works nicely for the tadpole expectation value as is evident in Fig. [5] The results at N^i = 2 
almost reproduces the analytical curve. So, the strategy to cure the pathology might seem promising. 

However, the gauge field expectation value scarcely benefits from this improvement procedure. We show the results 
in Fig. [5] for the a-fixed calculation (left) and the iJ^i-fixed calculation (right). The improvement seems rather better 
for the IR behavior with a fixed. The UV behavior with Ra fixed barely moves with increasing N^i- This is because 
the reference point for the i?^-fixed calculation is chosen at L = 700 where the deviation between the analytical and 
numerical results is acute as perceived from Fig. [3l So, must be comparable to ~ 700 to achieve a nice deal of 
improvement for the i?A-fixed results. We would observe better convergence if g^fiaL is smaller than 0.17 x 700 that 
we chose here. 



VII. DISCUSSIONS - FROM A "MODEL" TO A "THEORY" 



From the comparison between the analytical and numerical outputs, we can learn an important lesson; the numerical 
implementation like Eq. (jlip needs more and more caution as we approach the continuum limit in the transverse plane. 
This is a sort of dilemma. We should insert more and more sheets along the longitudinal direction as in Fig. 2] when 
we want to make use of a finer lattice or a larger volume in the simulation. Of course, the computation time increases 
as gets larger because the number of the integration variables is proportional to as seen in Eq. ([25)1 . 

We shall apply our results to the comparison of the initial energy density between the analytical formula by Eq. ([8|) 
and what was reported in Ref. (2^ . As we have discussed, the numerical formulation with the approximation ()lip 
underestimates the gauge field expectation value by a factor 0.546/0.15 = 3.64. Therefore, the initial energy density 
obtained in the numerical method is smaller than the analytical estimate by a factor 3.64^ — 13.2. 

We plot the initial energy density in Fig. [71 It should be mentioned that, though the important message from 
Ref. [22] is that the initial energy density at r = is ill-defined and logarithmically divergent in the Ra — > oo limit, 
our analytical calculation ends up with a finite value with the IR and UV cutoffs (i.e. a > and L < oo), which is 
just the same as the lattice discretized results. 
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FIG. 6; Gauge field expectation value for the SU(3) case at TV,, — 1 (same as shown previously), A'^,, = 2, and A'^,, = 10. The 
convergence to the analytical answer is much slower than the tadpole case especially for larger L. 
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FIG. 7: Initial energy density in the analytical calculation shown by the solid curve and the numerical estimate shown by the 
open circles. The cross points represent the data presented in Ref. 22] . 



It is notable that our numerical estimate is close to the data calculated in the numerical simulation in Ref. |22l | 
which is overlaid on Fig. [7] by the cross points. The small discrepancy should be explained by difference between the 
naive discretization in this work and the gauge-invariant lattice formulation in terms of the link variables. 

One should not jump to a conclusion, however, that the initial energy density becomes 13.2 times larger than the 
estimate in the previous numerical works. To earn a right interpretation, we have to understand how to specify the 
MV model parameter fi. The determination of fl has been carefully explained in the pioneering works [18. , ^0] and 
their choice {g — 2 and g^jl — 2 GeV for RHIC) has become a standard. Section V A in Ref. _2C| is actually devoted 
to elaborating that p, — 0.5 GeV can reproduce the total multiplicity at RHIC and then the following energy density 
is in a reasonable range. This procedure would make the discrepancy between the analytical and numerical results 
hidden behind the fitting; not only the energy density but also the multiplicity is accompanied by the factor 3.64. [The 
multiplicity is an integration over the gauge field expectation value divided by the particle dispersion relation.] That 
means that the value of ft would become •\/3.64 = 1.9 times larger in order to fit the total multiplicity without the 
overall factor 3.64, leading to 3.64^ = 13.2 times larger energy density. Consequently, even though the dimensionless 
coefficient has a substantial factor, the dimensional quantity in the physical unit remains unchanged. In other words 
the saturation model has only one dimensional scale fi and all the physical quantities are expected to scale with /2. 
Once p, is fixed by means of one of the experimental data set, other observables from the model calculation should be 
all consistent with the whole data set. 

Then, one might be confused at what the point of our finding is in this work. Our aim is, as emphasized in 
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FIG. 8: Schematic chart for the theory and model approaches. 



Introduction, to establish a correct theoretical procedure for the MV model. For that purpose the phenomenological 
success is inadequate. We must treat the MV model scale /I not as a fitting parameter but as an input deduced 
theoretically from a given Bjorken's x. This is, in principle, possible with the BFKL equation from which the saturation 
scale Qs is inferred as a function of x. [The nonlinear extension of the BFKL equation, namely the JIMWLK equation, 
is not necessary for the determination of Qs{x) as concisely explained in Ref. (26j.] The theoretical uncertainty by 
higher order corrections is still large, however, and the Golec-Biernat-Wiisthof fit multiplied by the atomic number 
factor A^/^ is usually regarded as a more trustworthy form of Qs{x). Then, within the framework of the MV model, 
the saturation scale can be defined in terms of ji as 



47r 



(29) 



where A is the IR cutoff. That is, A is given by Ra in the present work, while it should be the confinement scale 
~ 1 fm in reality. We can specify fl corresponding to relevant x through Qs{x) using Eq. 

We would call the above-mentioned strategy as the theory approach in contrast to the model approach where fi is 
phenomenologically fixed by the total multiplicity. Our point of view is summarized as a chart in Fig. [S] This present 
work is a first-step attempt to complete the full theory approach. 

Our finding that the longitudinal structure in infinitesimal extent has a significant effect may well be important 
also for the instability analyses in Ref. [1^, as we already mentioned in Introduction. Because the instability takes 
place with respect to the longitudinal fluctuations, it is naturally expected that the proper treatment of longitudinal 
randomness should alter the quantitative features discussed in Ref. [l^. It is an intriguing question whether the 
instability would remain weak as concluded in Ref. [25*1 or not with longitudinal randomness taken into account. If 
not, the instability might be fast and strong enough to account for the early thcrmalization. We have to emphasize 
that our iV^ is not the same as L^, employed in Ref. 25| in which the authors introduced the longitudinal coordinate 
to solve the equation of motion in three-dimensional space, but did not store the random sheets within infinitesimal 
extent in the initial condition. The vital difference lies in the fact that the number of the Monte-Carlo integration 
variables Pn{x±) {n — 1, . . . N,^) becomes greater with increasing 7V^. We conjecture that the Glasma instability would 
become stronger with iV„ > 1, but the quantitative analyses have to wait for at least times expensive computations 
as compared to Ref. [231, which is the problem to be investigated in the future. 
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